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The identification of pre-acceleration mechanisms for cosmic ray ions in supernova remnant shocks 
is an important problem in astrophysics. Recent particle-in-cell (PIC) shock simulations have shown 
that inclusion of the full electron kinetics yields non-time-stationary solutions, in contrast to previous 
hybrid (kinetic ions, fluid electrons) simulations. Here, by running a PIC code at high phase 
space resolution, ion acceleration mechanisms associated with the time dependence of a supercritical 
collisionless perpendicular shock are examined. In particular the components of J F-vdt are analysed 
along trajectories for ions that reach both high and low energies. Selection mechanisms for the ions 
that reach high energies are also examined. In contrast to quasi-stationary shock solutions, the 
suprathermal protons are selected from the background population on the basis of the time at 
which they arrive at the shock, and thus are generated in bursts. 

PACS numbers: 98.38.Mz 52.65.Rr 98.70.Sa 



I. INTRODUCTION 

Understanding the initial acceleration mechanisms for 
Galactic cosmic rays remains an outstanding problem in 
astrophysics. From energy balance considerations, super- 
nova remnants (SNRs) provide the most likely kinetic en- 
ergy source to sustain the cosmic ray population. The lo- 
cal acceleration of electrons has been indirectly observed 
at the expanding shock front of SNRs (see, for example, 
Ref. pj). However protons form the majority constituent 
of Galactic cosmic rays, and until recently observational 
evidence to link SNRs to local ion acceleration has been 
lacking. X-ray and 7-ray data from supernova remnant 
RX J1713. 7-3946 Q show energy spectra that can only 
be explained by accelerated ions. Several mechanisms are 
postulated to accelerate particles at SNR shocks. Fermi 
acceleration Q , which arises as a particle repeatedly scat- 
ters off turbulent structures on either side of the shock, is 
in principle capable of accelerating ions to high energies 
0. However, to work effectively an initial suprathermal 
population is required so that particles can re-cross the 
shock front ||. The identification and analysis of pre- 
acceleration mechanisms that can select and initiate the 
energisation of completely non-relativistic ions at SNR 
shocks from the background plasma is the subject of this 
paper. 

The Rankine-Hugoniot relations can be used to 
determine the discontinuity in bulk plasma parameters 
across a collisionless shock; that is, a shock where the 
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particle mean free path is much greater than length scales 
of interest. These relations are derived by applying the 
magnetohydrodynamic (MHD) conservation equations in 
the far upstream and downstream limits, away from the 
shock. Further conditions are imposed by the fact that a 
shock must also increase entropy, so that no subsonic flow 
can spontaneously become supersonic. For an Alfvenic 
Mach number Ma > 3, the shock is supercritical in that 
the increase in entropy, and in ion heating, required by 
the Rankine-Hugoniot relations is achieved via the ion 
kinetics; at least in part, by reflection of a fraction of 
upstream ions at the shock. The generic supercritical, 
quasi-perpendicular, collisionless shock in which ions re- 
flect and gyrate in a foot region upstream has been sug- 
gested by hybrid (particle ions, fluid electrons) simula- 
tions (see, for example, Refs. |7j,|8|,|9j), and confirmed by 
in-situ observations of the Earth's bow shock [Toj . 

To study the acceleration of ions and electrons, a 
fully kinetic treatment is in principle necessary for both 
species, and this can be closely approximated by particle- 
in-cell (PIC) techniques. Physical mechanisms operating 
on electron kinetic lengthscales and timescales are sig- 
nificant both for aspects of macroscopic structure (for 
example, the shock ramp width scales as c/u> pe ), and for 
microscopic processes affecting the ions (such as caviton 
formation and dissolution) . Whether such effects are im- 
portant in any given scenario can be estimated, to some 
extent, by analytical means, as we discuss in detail be- 
low. Importantly, however, it is known that inclusion of 
the full electron kinetics can significantly alter the dy- 
namics of the shock. For example, hybrid simulations for 
certain parameters 0, El produce time-stationary shock 
solutions, whereas for the same parameters PIC simula- 
tions reveal a dynamic, reforming, shock structure. Fur- 
thermore the extent to which an individual ion responds 
to phenomena on electron kinetic scales must depend on 
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that ion's cyclotron radius, and hence its energy. It fol- 
lows that for studies of ion acceleration at shocks, as 
in the present paper, retention of full electron kinetics 
is desirable in order to resolve fully the shock dynam- 
ics (see also Refs. 0, an d the ion dynamics. We 
have previously presented results of PIC code simula- 
tions |l2l fl3L that have high resolution in real space 
and phase space, over relatively long run times, for pa- 
rameters relevant to shocks at supernova remnants. Our 
most recent results 01 show that the time-dependent 
electromagnetic fields at the reforming shock can accel- 
erate inflow ions from background to suprathermal ener- 
gies. This provides a source population which may sub- 
sequently be accelerated to produce high energy cosmic 
rays. 

In the present paper, we focus on the specific nature of 
this ion acceleration mechanism. This requires careful ex- 
amination of the physics of the interaction between parti- 
cles and fields as they evolve over time. We first introduce 
a methodology for simplifying the raw data, obtained 
on spatio-temporal scales spanning those of the electrons 
and ions, into data suitable for examining events on the 
spatio-temporal scales of interest for ion acceleration. We 
then examine the detailed dynamics of the ion interac- 
tions with the shock front, including a comparison be- 
tween ions that eventually reach the highest and lowest 
energies downstream. The time at which particles are 
incident on the temporally evolving shock structure is 
found to be a key discriminant in the subsequent energi- 
sation. 



II. SIMULATION METHOD 

The technical basis of the simulations was recently re- 
ported in Ref. 0; let us reiterate briefly for complete- 
ness. We use a relativistic electromagnetic particle-in-cell 
(PIC) code to simulate the structure and evolution of 
a supercritical, collisionless, perpendicular magnetosonic 
shock. In a PIC simulation the distribution functions 
of all particle species are represented by computational 
super-particles, whilst the electromagnetic fields are de- 
fined on a spatial grid. Particle trajectories are evolved 
from the fields via the Lorentz force, then the fields are 
evolved from the new particle positions via Maxwell's 
equations . The code we use to simulate the shock is 
based on that described in Ref.[l^|. and has been used 
recently to examine electron and ion acceleration in SNR 
shocks [HEl Q E3 All vector fields, bulk plasma 
properties and particle velocities arc functions of one spa- 
tial co-ordinate (x), and time. This simplification en- 
ables detailed phase space resolution for relatively long 
run times, however it constrains magnetic fields: since 
V • B = we have B x = constant, taken as zero here 
in strict perpendicular geometry. PIC simulations in two 
spatial dimensions (see, for example, Ref. ^^|) that re- 
lax this constraint yield overall shock dynamics that are 
consistent with the results seen here. 



We present results from simulations of a perpendic- 
ular shock propagating into a magnetic field (B z ^\) of 
1 x 10~ 7 Tesla, a value consistent with those expected at 
supernova remnants |l9j . The ratio of electron plasma 
frequency to electron cyclotron frequency uj pe /uj ce = 20, 
and the upstream ratio of plasma thermal pressure to 
magnetic field pressure, (3 = 0.15. The shock has an 
Alfvenic Mach number (Ma) of 10.5, and the simulation 
mass ratio for ions and electrons Mr = mi/m e = 20, in 
common with Refs.m* El> El HiJ- This reduced mass 
ratio is necessary to enable ion and electron time scales 
to be captured within the same simulation, with reason- 
able computational overheads. Previous PIC simulations 
for physical, and a range of non-physical, uii/m e show a 
variety of kinetic instabilities in the foot region [ill l2lj . 
Here we find that the ion dynamics are insensitive to 
structures on electron scales, associated with these insta- 
bilities. 

We also follow Refs. [II, E3, H [lij in using the piston 
method (see, for example, Ref. [ij, and references therein) 
to set up the shock. Particles are injected on the left 
hand side of the simulation box with a drift speed Vi n j , 
modified by a small random velocity drawn from a ther- 
mal distribution, characterised by u t herm- At the particle 
injection boundary, the magnetic field (B z a) is constant 
and the electric field (E y ^) is calculated self consistently. 
The right hand boundary is taken to be a perfectly con- 
ducting, perfectly reflecting wall. Particles reflect off this 
boundary, and a shock then forms, propagating to the left 
through the simulation box; sufficient time is allowed for 
the shock to propagate sufficiently far upstream that the 
boundary conditions arc no longer important. The size 
of a grid cell is defined as a Debye length (Xd), and the 
time ste p is set to less than A^> / c, for numerical stability 
reasons [23. To enable the shock and particle dynamics 
to be followed over extended time-scales, whilst retaining 
high particle density, a simple shock following algorithm 
is implemented. This holds the peak in magnetic field at 
8A ci from the left-hand boundary (for details see Ref. 0] , 
Appendix A). This distance is chosen so that no particles 
that are reflected off the shock subsequently reach the up- 
stream boundary, whilst it leaves a significant region of 
the simulation box (around 20A C i) downstream. 



III. RESULTS 

Full simulation of the non-time-stationary features of 
a collisionless shock requires the retention of electron dy- 
namics; see, for example Ref. @- However, resolving 
features on electron scales also introduces processes that 
do not couple strongly to the processes that operate on 
ion scales, which are the focus of the present paper. For 
example, the observed electron scale electrostatic soli- 
tary wave structures can lead to electron acceleration 
[l4| . but do not significantly affect the ions. As an aid 
to interpreting the interactions occurring within a PIC 
simulation that give rise to ion acceleration, we present 
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a framework which distinguishes structure and dynam- 
ics on electron kinetic scales from those relevant to ions. 
The ion trajectories that we present here are however all 
evolved self-consistently within the full electro-magnetic 
fields of the PIC simulation 



Electric Potential on Ion Spatio- Temporal 
Scales 



Resolving the full electron and ion kinetics in the PIC 
simulation establishes two distinct spatio-temporal scales 
on which physical processes can occur. On sufficiently 
fast time scales and short length scales there are, for 
example, plasma oscillations that lead to fluctuations 
in charge density. However, on longer spatio-temporal 
scales the plasma is quasi-neutral but still supports bulk 
electric fields. Wc can obtain an expression for these bulk 
fields by treating the plasma as two fluids, ions and elec- 
trons (for a more general multi-fluid treatment see, for 
example Rcfs. |23|, |2 113 ) > governed by the momentum 
equations 



: Dt 
1 Dt 



-en e (E + v e AB) - VP e - n e v el y e , (1) 
q i n i (E + v i AB)-VP i -n i u ie v i . (2) 



The final terms in Eqs. Q and represent momentum 
transfer between species via forces not included in the 
macroscopic fields. Here, we can assume that on the 
time-scale of ion interaction with the shock, these terms 
are negligible for the ions. 

We wish to consider space and time varying electro- 
magnetic fields that only affect the ions, so that on 
ion scales we can take the limit in which the elec- 
trons respond instantaneously as a charge neutralising 
fluid. This implies a vanishing electron inertial term, the 
"massless electron fluid" limit often used in hybrid codes: 
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We neglect the electron pressure gradient because it is 
significant on electron, rather than ion scales, however we 
anticipate that this approximation will be least reliable 
in the shock ramp. 

We can relate v e directly to v.; via the current. On 
the spatio-temporal scales on which the electron-proton 
(?« = e) plasma is quasi-neutral (n% ~ n e = n), 

J ~ en(v, - v e ). (4) 

Substitution for v e from Eq. into Eq. (|TJ then gives 

J 



~ E 



AB. 



(5) 



Consistent with this low frequency approximation we ne- 
glect the displacement current (the non-radiative limit), 



giving Ampere's law 

V A B = /i J. 
This implies the standard expression 
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Together, Eqs. © and Q give 
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which can then be rearranged to yield E. 

In the one dimensional geometry of our simulation 
Eq. (JSJ can be simplified by noting V = (c? x ,0, 0), 
thus (B • V)B = 0. Further simplification arises if wc 
note that generally in our simulations, v z << v y , thus 



(vi, y B z - v hZ B y ) 
Eq. JSJ then gives 
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B z . Rearranging a simplified 
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Here E x ,i is the x component, and E y j the y component, 
of the electric field on the slow, ion spatio-temporal scales 
on which the plasma is quasi-neutral. 

Substitution of our simplified, ion scale, electric field 
from Eq. JHJ), into the ion force equation Eq. @, leads to 
an expression whose x component is 
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It follows from Eq. that the bulk force on the ion 
fluid is due to gradients in magnetic and plasma pressure. 
The potentials (which act on individual particles) follow 
from Eqs. @ and ifTUj). 

Figure ^ demonstrates the extent to which this ap- 
proximate analytical treatment provides a guide to the 
ion behaviour that is calculated from first principles in 
the PIC code. It compares the time evolution of the 
potential, <f> = J E x dx, obtained directly from the PIC 
code, to that calculated on ion scales, 4>i = / E Xt idx, us- 
ing Eq. ©. The path chosen for the spatial integration 
is that of an ion that reaches a high energy on leaving 
the shock front. Fig. ^ demonstrates that E x ^ defined in 
Eq. I© is a useful guide, and hence the analysis above 
captures much of the key physics. The ion scale bulk 
potential essentially averages over the small scale fluctu- 
ations of the "raw" potential. We can see from Fig. ^ 
however, that the average values of 4>i depart from that 
of the full potential </> where the ion interacts with the 
shock ramp: first during a reflection at t = 3.5 — 3.7, 
and during a subsequent transmission to downstream at 
t = 5 — 5.2 x 2ituj~ , In the discussion below, we will 
calculate the ion energetics from the full electromagnetic 
fields of the PIC simulation. 
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B. Ion Acceleration 

To study the physical processes that cause ion acceler- 
ation, we evaluate the changes in kinetic energy of ions 
during their interaction with the shock. Here the ion 
Lorentz factor 7 ~ 1, therefore we can neglect relativis- 
tic effects. We have 

F = m$=g(E + vAB) (12) 
dt 

where, in our collisionless plasma, E and B in the Lorentz 
force refer to the fields in the PIC simulation. Thus 

F ' v= IG mv2 ) =9E - v (i3) 

Integration along a computed ion trajectory then implies 
that the kinetic energy acquired is: 

-mv 2 =q E • vdt. (14) 

^ J trajectory 

1. Highly energetic ions 

Previous PIC simulations 01 have shown that the 
downstream proton population has a continuous distri- 
bution of energies from zero up to ~ 6 times the ion 
injection energy, £i n j = ^mvf n j, in the frame in which 
the downstream plasma is at rest. We now examine the 
dynamics of these ions in more detail. Figure El presents 
the results of evaluating Eq. ((HI) for a selected group 
of protons that become highly energised. The top panel 
(panel 1) displays the kinetic energy over time, calcu- 
lated from J qE ■ vdt along the particle trajectory. Panel 
2 shows only the ^-component, J qE x v x dt normal to the 
shock, and panel 3 shows only the y-componcnt along the 
shock front. The z-component is omitted as it remains 
identically zero, due to the configuration of the simula- 
tion domain. Panel 4 displays the potential, <fi = J E x dx, 
computed directly from the PIC code, at the x-position 
of the ions at the current time. Panel 5 shows the y- 
positions. In the lowest panel (panel 6) the x-positions 
are shown in relation to the spatio-tcmporally evolving 
potential structure on ion scales, <pi = J E x ^dx computed 
usin g Eq . @. Comparison of this panel with Fig. 8 of 
Ref. [12] shows that 4>i captures the qualitative features 
of the electromagnetic fields. To complement this infor- 
mation, Fig. |21 shows the trajectory of a high energy ion 
in the x-y plane (in the simulations we evolve the three 
components of the particle velocity v(x, t), and these can 
be integrated to provide a trajectory in three dimensional 
configuration space). As with all results in this paper, 
data is presented in the downstream rest frame, and has 
been obtained from a segment of the simulation when 
the shock is propagating independently of the boundary 
conditions. Units are normalised to the upstream ion pa- 
rameters, that is, X c i the upstream ion cyclotron radius, 
and u>d the upstream ion cyclotron frequency. 



Figure [3 shows that the ions that become highly ener- 
gised remain close in phase space throughout their inter- 
actions with the shock. After passing through the shock, 
local fluctuations in the fields lead to some divergence in 
the y-component of J q~E-vdt and the y-position. Panel 6 
shows the shock propagating in the negative x-direction 
over time, while undergoing reformation cycles. The size 
and depth of the potential well varies over the course of 
a reformation cycle, on a time scale comparable to the 
local ion cyclotron period, as discussed in detail by Lee 
et al. 

If we follow the path through the shock region of an in- 
dividual ion that eventually reaches high energy, a series 
of events occurs. The ion is initially co-moving with the 
plasma at the inflow speed. This corresponds to the lin- 
ear increase in x-position (panel 6), with no translation 
in the y-direction (panel 5 and Fig- EJ) , all at the inflow 
energy (panel 1). It should be noted that in panel 1, the 
energies are initially zero, because the start of the inte- 
gration path for J q~E ■ vdt is chosen to be just inside the 
simulation domain, where the ions are already co-moving 
with the plasma. This location is upstream of the shock 
foot, and the potential in this region is therefore constant, 
with the only variation due to high frequency fluctuations 
in E. 

After approximately 3.25 ion cyclotron periods, the ion 
enters the potential well upstream of the shock front, 
point a in Fig. |21 At this time, the shock is most fully 
formed, the shock jump is close to maximal, so that the 
dB 2 /dx term in the potential from Eq. © is close to 
maximal also. The x-component of J q~E-vdt shows a de- 
crease in energy as the ion journeys further into the well 
(panel 2), with a corresponding decrease in kinetic en- 
ergy (panel 1); this follows from a negative value for E x . 
Between points a and b there a is decrease in x-vclocity 
as the ion gets closer to the shock and the magnetic field 
increases, which is accompanied by a weak drift in — y, 
see also Fig. |3 

At the time corresponding to point 6, the kinetic en- 
ergy reaches a minimum (panel 1), when the ion is near 
the deepest point in the potential well (panels 4 and 6). 
The ion has now stopped its progression toward the shock 
front (Fig. |3J), and reflection back into the upstream re- 
gion has begun (panel 6). 

By time c the ion has begun to climb back out of the 
potential well, away from the shock front (panel 4). After 
c, drift in +y begins as the ions move in —x into the foot 
region, see also Fig. |21 Meanwhile the x-component of 
J q~E ■ vdt starts to increase rapidly (panel 2), due to the 
strength of E x (shown by the gradient of (f>, panel 4). 

By point d the particle has moved back to the extreme 
upstream edge of the potential well. It then remains 
moving along a contour of <fi ~ (panels 4 and 6) for a 
time approximately equal to one upstream ion cyclotron 
period (27ro;~ 1 ). Between d and e the value of E x local 
to the ion is lower, so that the associated energisation 
rate is also less; however the positive y-componcnt of the 
gyromotion continues (panel 5 and Fig.|3J), and since the 
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motional E y is positive, this gives an energy gain in the 
y-component A£ y = Amv^/2 (panel 3). The gyromotion 
of the particle eventually leads to a positive ^-component 
of velocity (Fig. EJ , so that at point e the particle finally 
leaves the extreme upstream edge of the well and passes 
through the potential well for a second time (panels 4 
and 6). This marks the end of a prolonged episode of 
energy gain, which now stops as y-drift ceases (panels 1 
and 3), and the particle settles into its stable downstream 
gyromotion (Fig. EJ. The ions cross the saddle in the 
potential 4>(x, t) (as shown in panel 6), leading to a brief 
energy loss then gain via the x-component of J qE • vrfi 
(panel 2), before the ion passes thorough the shock front 
(point /, panels 3 and 6), and gyrates away from the shock 
into the downstream region. The ion energy now exceeds 
its initial value by a factor of approximately six. 

In summary there arc two stages of acceleration as 
shown in Fig. EJ normal reflection from the temporar- 
ily stationary shock front into the foot region, followed 
by energisation during transverse drift across the shock 
front. 



2. Weakly energised ions 

Further insight into the energisation process can be 
gained by comparing trajectories for ions that become 
highly energised (Fig. EJ), to those for a group of ions 
that have only low energies (< £i n j) on finally entering 
the downstream region fFig. Hjl . and remain in the core 
of the downstream particle distribution. 

The trajectories for the ions that are not subsequently 
energised are initially similar to those for the ions that 
eventually reach the higher energies. The primary dif- 
ference is the timing of their first encounter with the 
shock. We have identified two distinct groups of low en- 
ergy ions, and these arc shown in Fig. El The first group 
arrives at the shock front just as the shock is advancing 
(t = 3.1 x 2ttuj~ 1 ), and the second when the shock is 
decaying (t = 4x 27rw t J i 1 ). 

The reforming shock progresses upstream (downwards 
in panel 6 of Figs. and 0} in a stepwise fashion. The 
first group of ions (upper trajectories in panel 6 of Fig.EJ 
encounter the advancing shock when the shock jump is 
sufficiently large to cause reflection (between points A 
and B). Their trajectories up to this point are akin to 
the trajectories up to point b in Fig. 0for the ions that 
have become highly energised: they have entered the foot 
region (point A) and been deflected in the — y direction, 
whilst losing energy via a decrease in J qE x v x dt. How- 
ever, at this time the shock speed is maximal, so their 
velocity component in —x is smaller than that of the 
shock itself, and the shock overtakes them. They then 
co-move with the shock front for about an upstream ion 
cyclotron period, before moving downstream. 

On the other hand the second group of ions encounter 
the shock when the potential is decaying (point A'). 
They then pass through the potential well (point B'), 



and reach the shock without reflection, where their v±_ in- 
creases (Fig. EJ) along with the bulk B field downstream. 

Regardless of their energisation or of the details of 
their dynamics, the guiding center velocity of all ions 
goes to zero once they have propagated sufficiently far 
downstream of the shock front. This is to be expected, 
because the far downstream frame defines the rest frame 
of the plasma, as noted previously (Figs. El and EJ)- 

Of further interest is the y-motion of the two groups of 
low energy ions (panel 5 of Fig.^J. Those that enter the 
shock before the high energy ions have little movement 
of their gyrocenters in y, but those that enter after the 
high energy ions, have a significant — y drift velocity; see 
also Fig. Both these patterns are in contrast to the 
high energy ions that have gyrocenters that drift in the 
+y direction (panel 5 of Fig.0 and Fig. EJ). 

Finally we can compare the kinetic energy gain com- 
puted directly from the electromagnetic fields of the PIC 
simulation, with that given by the fields on ion scales 
estimated from Eqs. (JHJ) and llOjl . In Fig. El we plot the 
total kinetic energy as a function of time along the trajec- 
tory of an ion that is reflected and reaches suprathermal 
energy. This plot shows a close correspondence between 
the two curves, suggesting that the fine structure on elec- 
tron scales does not affect the final energy gain of these 
ions. However, as we have seen from Fig. there is a 
discrepancy between the x-componcnt of the electric field 
at the shock ramp obtained from the simulation directly, 
and from Eq. (JHJ) . It therefore follows from Fig. El that 
the value of the shock ramp does not strongly affect the 
overall energy gain of these ions. This is consistent with 
energisation being associated with electromagnetic fields 
away from the ramp, the role of the ramp potential be- 
ing simply to reflect the ions. However, details of the 
energetics of low energy ions that are not reflected, may 
depend on the value of the shock ramp potential. 



C. Role of Shock Reformation 

Having examined the trajectories of ions that reach 
a variety of energies, let us now examine the selection 
mechanisms that give rise to different histories and en er- 
gisation. For a time stationary shock, iBurgess et alJ Q 
examined the origins in phase space of ions that even- 
tually reach differing energies. Particles from the ex- 
terna of the velocity space distribution upstream of the 
shock were found to be preferentially reflected further 
upstream, and so energised to higher energies; whereas 
ions from the core of the distribution passed through the 
shock, moving little or no distance upstream. To estab- 
lish whether the same selection mechanism is at work 
in our dynamic reforming shock, we have constructed in 
Fig. a series of plots that may be compared with Fig. 1 
of Ref. Q . Figure shows the ion phase space (v x and 
v y vs. x) for groups of ions at differing initial perpen- 
dicular velocities, at a time when the potential well is at 
its narrowest, and the reformation cycle has just com- 
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menced, corresponding to t — 4 x inujZ^ 1 in Fig. [2J In 
contrast to the results obtained in Ref. [!j for the case of a 
time stationary shock, we find that at a reforming shock, 
ions that are initially in the core of the distribution, as 
well as those from the tails, are reflected back into the 
foot region. Also, the distance travelled back into the 
foot region (and hence the energy gained) appears inde- 
pendent of the initial perpendicular velocity of the ions. 
Whether or not a given ion is reflected depends on its 
normal velocity in comparison with the time-dependent 
shock potential. Thus the timing at which ions arrive at 
the shock front determines their final location in velocity 
space. 

Overall, examination of the shock dynamics in relation 
to ion trajectories shows that the ions that are ultimately 
highly energised (Fig. [2J) reflect from the shock front just 
as it becomes stationary, and pass through the foot region 
saddle in (f>(x,t) (panel 6). The ions that meet the shock 
front prior to this (upper group in Fig. 0J) interact with 
a shock front that is moving rapidly forward through the 
simulation domain, and so do not gain sufficient velocity 
to outpace it. The ions that interact later (lower group 
in Fig. meet a weakening shock front with a wider 
potential well, so that they are not reflected. In contrast 
to Fig. the ions in Fig. 0] experience neither an x- 
energy gain on moving back to the upstream side of the 
potential well, nor a sustained period of y-energy gain 
as they subsequently co-move with the upstream edge of 
the well. 

The present results suggest that the time-evolving 
shock dynamics, and in particular the timing of the in- 
teraction between ion and shock, govern the selection 
process determining which ions undergo pre-acceleration 
into a suprathermal population that may subsequently 
become cosmic rays. 

IV. CONCLUSIONS 

We have examined in detail the dynamics of suprather- 
mal ions generated in PIC code simulations of quasi- 
perpendicular reforming shocks. Importantly, this en- 
ergisation is not found in stationary shock solutions. We 
find that: 

1. The shock structure reforms on a timescales of the 
order of the local ion cyclotron period. This is 
shown clearly if the electromagnetic fields are cast 
in the form of a potential, after removing small 
scale effects, to leave only terms relevant on ion 
spatio-temporal scales. 

2. The time-dependence of the shock dominates the 
selection of which ions are accelerated to suprather- 
mal energies. Ions that reach the shock when its 
ramp, and hence potential, are maximal, are re- 
flected and subsequently gain energy by drifting in 
the time-dependent fields tangential to the shock 
front. 



3. This selection is in contrast to a time-stationary 
shock, where the selection mechanism depends 
upon the initial ion velocity perpendicular to the 
magnetic field, those ions coming from the tails of 
the distribution being preferentially reflected, and 
so energised. 

4. These factors lead to high energy ion creation oc- 
curring in bursts. 

The present simulations are conducted in a lx3v geom- 
etry at a low m,i/m e . Whilst the work in Ref. 0| shows 
no major differences in higher spatial dimensions for PIC 
simulations, hybrid simulations in 2x3v, for longer run 
times, show that the current that can now exist alongthe 
shock front can lead to current-driven instabilities [26j. 
These instabilities may act to change the shock structure 
in the y-direction, and so alter ion and electron dynamics 
across the shock front, affecting both ion and electron ac- 
celeration. Simulations with more realistic mi/m e ratios 
[ill l2lj show alterations in the electron scale physics in 
the foot region. However, we have shown that the elec- 
tron scale physics has little effect on ion spatio-temporal 
scales. 

The plasmas simulated here are pure hydrogen, in that 
there are only two species, protons and electrons. Both 
species have a Maxwellian distribution with the same 
temperature. The addition of pickup ions to the sim- 
ulation, for example in relation to the heliospheric ter- 
mination shock, where hybrid simulation have already 
been carried out [53, HU, would allow the acceleration 
processes relevant to anomalous cosmic ray production 
to be examined. This will be the subject of future work. 

The fundamental plasma physics processes underlying 
the ion acceleration from background to suprathermal 
energies (10 to 20 MeV), reported in the SNR shock sim- 
ulations of Ref. 0| , have been elucidated in the present 
paper. Specifically, we have explained the role of, and 
interplay between, the key elements anticipated at the 
end of the Appendix to Ref. 0|. We have shown that, 
while an electron fluid approximation captures some of 
the key physics, the shock reformation dynamics arising 
from our fully kinetic PIC treatment are central to the ion 
acceleration mechanism. This work provides a clear first 
principles explanation for the ion acceleration that is ob- 
served in our simulations, which appears to be a strong 
candidate injection mechanism for Galactic cosmic ray 
protons. 
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FIG. I: <f> = f E x dx from PIC code (gray), and 4>i = f E x ^dx 
calculated from Eq. @ (black) along the trajectory of a high 
energy ion. In calculating E Xti we compute the ion flow ve- 
locity from the mean velocity of all ions within 0.02A C j ~ 
A C e/2 ~ 10 grid cells of the particle position. 



FIG. 2: Trajectories for 4 ions that reach high energies. Panel 
1 shows f qE ■ wdt along each trajectory. Panel 2 the x- 
component of f qE ■ vdt, panel 3 the y-component. The z- 
component is not shown as it remain identically 0. Panel 4 
displays <j> = J E x dx. Panel 5, y-position and panel 6, re- 
position plotted over the potential on ion scales, derived from 
Eq. here the shock is propagating towards lower values 
of x as t increases. The vertical lines correspond to times of 
change during the black trajectory. 



FIG. 3: Position x vs. y for a high energy ion. This ion 
follows the black trajectory in Fig. [5] and the timing points 
a-f are indicated. 



FIG. 4: Trajectories for four ions that remain at low energies 
when crossing the shock front. Panels and colors are as Fig. 01 



FIG. 5: Position x vs. y for two low energy ions. The black 
and gray trajectories from Fig. 21 are represented here, with 
the timing points A, B, A' , B' indicated. 



FIG. 6: Kinetic energy calculated from Eq. 11141 using the PIC 
simulation E (gray) and that on ion scales Ej (black) from 
Eqs. © and JTDJ, for the black ion in Figs. and El The 
bulk ion velocity Vix,y in Eqs. © and 1101 is calculated from 
the mean velocity of the ions within 10 grid cells ~ 0.02A c i ~ 
A ce /2 of the particle position. 



FIG. 7: Phase space plots of v m (left) and v y (right) vs. x 
for groups of particles from differing regions of initial velocity 
space, at the instant t = 4x 27TLJ~ t . 
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FIG. 1: <f) = J E x dx from PIC code (gray), and 4>i = J E x ^dx 
calculated from Eq. @ (black) along the trajectory of a high 
energy ion. In calculating E m> i we compute the ion flow ve- 
locity from the mean velocity of all ions within 0.02A C ; ~ 
A ee /2 ~ 10 grid cells of the particle position. 
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FIG. 2: Trajectories for 4 ions that reach high energies. Panel 
1 shows J gE ■ wdt along each trajectory. Panel 2 the x- 
component of f qE ■ wdt, panel 3 the y-component. The z- 
component is not shown as it remain identically 0. Panel 4 
displays <f> = j E x dx. Panel 5, y-position and panel 6, re- 
position plotted over the potential on ion scales, derived from 
Eq. here the shock is propagating towards lower values 
of x as t increases. The vertical lines correspond to times of 
change during the black trajectory. 
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FIG. 3: Position x vs. y for a high energy ion. This ion 
follows the black trajectory in Fig. and the timing points 
a-f are indicated. 
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FIG. 6: Kinetic energy calculated from Eq. Ijl4^ using the PIC 
simulation E (gray) and that on ion scales Ej (black) from 
Eqs. © and 0, for the black ion in Figs. and 01 The 
bulk ion velocity Vi m>y in Eqs. JJJJ and llOU is calculated from 
the mean velocity of the ions within 10 grid cells ~ 0.02A c i ~ 
X ce /2 of the particle position. 
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FIG. 7: Phase space plots of v x (left) and v y (right) vs. x 
for groups of particles from differing regions of initial velocity 
space, at the instant t — Ax 2ivu)^f. 



